/*
 *  hyperelastic.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 1/04/11.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "hyperelasticstretches.h"

HyperelasticStretches::HyperelasticStretches()
:	Material("",""),
	lambda1(1.0),
	lambda2(1.0),
	lambda3(1.0),
	W1(0.0),
	W2(0.0),
	W3(0.0),
	W11(0.0),
	W12(0.0),
	W13(0.0),
	W22(0.0),
	W23(0.0),
	W33(0.0)
{}

// ***

HyperelasticStretches::HyperelasticStretches(const string& na, const string & t)
:	Material(na,t),
	lambda1(1.0),
	lambda2(1.0),
	lambda3(1.0),
	W1(0.0),
	W2(0.0),
	W3(0.0),
	W11(0.0),
	W12(0.0),
	W13(0.0),
	W22(0.0),
	W23(0.0),
	W33(0.0)
{}

// ***

void HyperelasticStretches::calcStretches( const VectorXd & C )
{
	double tol = 1e-12;
	double lam1, lam2, lam3;

	lam1 = C(0);
	lam2 = C(1);
	lam3 = C(2);

	if(fabs(lam1) <= tol && fabs(lam2) <= tol && fabs(lam3) <= tol)  //If matrix C is diagonal it's not neccesary to calculate the principals stretches.
	{
		lambda1 = sqrt(lam1);
		lambda2 = sqrt(lam2);
		lambda3 = sqrt(lam3);
		
		J = lambda1 * lambda2 * lambda3;
		    
		N1 << 1.0,
		      0.0,
		      0.0;
		
		N2 << 0.0,
		      1.0,
		      0.0;
		
		N3 << 0.0,
		      0.0,
		      1.0;
		
	}
	else  // Obtaining the principals stretches and principal directions
	{
		InverseVoigtFormat(CSym,C);
		
		SelfAdjointEigenSolver<Matrix3d> eigensolver(CSym);
		
		stretches = eigensolver.eigenvalues();
		ppaldirecctions = eigensolver.eigenvectors();
		
		lam1 = stretches(0);
		lam2 = stretches(1);
		lam3 = stretches(2);

		lambda1 = sqrt(lam1);
		lambda2 = sqrt(lam2);
		lambda3 = sqrt(lam3);
		
		J = lambda1 * lambda2 * lambda3;
		
		N1 = ppaldirecctions.col(0);
		N2 = ppaldirecctions.col(1);
		N3 = ppaldirecctions.col(2);
		
// 		cout << "N= " << endl << ppaldirecctions << endl;
	}

	if(!hypertype.compare("decoupled"))
	{
		Jm =  pow(J, - 1.0 / 3.0);
		
		lambda1m = Jm * lambda1;
		lambda2m = Jm * lambda2;
		lambda3m = Jm * lambda3;
	}
}

// ***

void HyperelasticStretches::calcInvariants( const VectorXd & C )
{
    calcStretches(C);

    I = lambda1*lambda1 + lambda2*lambda2 + lambda3*lambda3;
    II = pow(lambda1*lambda2,2) + pow(lambda1*lambda3,2) + pow(lambda2*lambda3, 2);
    III = pow(lambda1*lambda2*lambda3,2);
    
}

// ***

void HyperelasticStretches::calcCiTensor( const VectorXd & C)
{
    InverseVoigtFormat(CSym,C);

	CiSym = CSym.inverse();

	CiTensor(1,1) = CiSym(0,0);
	CiTensor(2,2) = CiSym(1,1);
	CiTensor(3,3) = CiSym(2,2);
	CiTensor(2,1) = CiSym(0,1);
	CiTensor(3,1) = CiSym(0,2);
	CiTensor(3,2) = CiSym(1,2);
}

// ***

void HyperelasticStretches::calcCTensor (const VectorXd & C)
{
	InverseVoigtFormat(CTensor,C);
}

// ***

void HyperelasticStretches::calcP(const VectorXd & C)
{
	InverseVoigtFormat(CTensor,C);
    calcCiTensor(C);
    P = Id4 - 1.0 / 3.0 * dyadic(CiTensor, CTensor);
}

// ***

void HyperelasticStretches::calcPt(const VectorXd & C)
{
	InverseVoigtFormat(CTensor,C);
    calcCiTensor(C);
    Pt = Id4 - 1.0 / 3.0 * dyadic(CTensor, CiTensor);
}

// ***

void HyperelasticStretches::calctildeP(const VectorXd & C)
{
    calcCiTensor(C);
    tildeP = selfOdot(CiTensor) - 1.0/3.0 * selfDyadic(CiTensor);
}

// ***

void HyperelasticStretches::calcW( const VectorXd & C )
{
		calcStretches(C);
		strainEnergy();
}

// ***

void HyperelasticStretches::calcS( const VectorXd & C)
{
	calcStretches(C);
	DstrainEnergy();
	
	if(!hypertype.compare("decoupled"))
	{
		calcSiso(C);
		calcSvol(C);
		S = Svol + Siso;
	}
	else
	{
		//   S = 2.0 * Sum_a(W_a / lambda_a * N_a * N_a)
		
		if(lambda2 == 0.0 and lambda3 == 0.0)	//means elasticlinear or elasticlog
		{
			S1 = W1 / lambda1;
			
			SMat = S1 * N1 * N1.transpose();
		}
		else
		{
			S1 = W1 / lambda1;
			S2 = W2 / lambda2;
			S3 = W3 / lambda3;
			
			SMat = S1 * N1 * N1.transpose() +  S2 * N2 * N2.transpose() + S3 * N3 * N3.transpose();
		}
		VoigtFormat(SMat, S);
	} 
}

// ***

double  HyperelasticStretches::getDDW(const VectorXd & C )
{
	int i, j;
	calcStretches(C);
	DDstrainEnergy();
	
	return W11;
}

//***

void HyperelasticStretches::calcSvol(const VectorXd & C)
{
	calcStretches(C);
	DstrainEnergy();
	
	SvolMat =  J * dWvol * ( N1 * N1.transpose() /  (lambda1 * lambda1) + N2 * N2.transpose() / (lambda2 * lambda2) + N3 * N3.transpose() / (lambda3 * lambda3) );
	
	VoigtFormat(SvolMat,Svol);
	
	//cout << "Svol= " << endl << Svol << endl;
}

//***

void HyperelasticStretches::calcSiso( const VectorXd & C)
{
	calcStretches(C);
	DstrainEnergy();
	
	Siso1 = ( 2.0 / 3.0 * lambda1m * Wiso1 - 1.0 / 3.0 * (lambda2m * Wiso2 + lambda3m * Wiso3) ) / (lambda1 * lambda1);
	Siso2 = ( 2.0 / 3.0 * lambda2m * Wiso2 - 1.0 / 3.0 * (lambda1m * Wiso1 + lambda3m * Wiso3) ) / (lambda2 * lambda2);
	Siso3 = ( 2.0 / 3.0 * lambda3m * Wiso3 - 1.0 / 3.0 * (lambda1m * Wiso1 + lambda2m * Wiso2) ) / (lambda3 * lambda3);
	
	SisoMat = Siso1 * N1 * N1.transpose() + Siso2 * N2 * N2.transpose() + Siso3 * N3 * N3.transpose();
	
	VoigtFormat(SisoMat,Siso);
	
	//cout << "Siso= " << endl << Siso << endl;
}

//***

void HyperelasticStretches::calcCt( const VectorXd & C)
{
	double tol = 1e-8;
	
	if(N1N1.size()==0)	N1N1.setZero(6);
	if(N2N2.size()==0)	N2N2.setZero(6);
	if(N3N3.size()==0)	N3N3.setZero(6);
	if(N1N2.size()==0)	N1N2.setZero(6);
	if(N1N3.size()==0)	N1N3.setZero(6);
	if(N2N1.size()==0)	N2N1.setZero(6);
	if(N2N3.size()==0)	N2N3.setZero(6);
	if(N3N1.size()==0)	N3N1.setZero(6);
	if(N3N2.size()==0)	N3N2.setZero(6);
	
	calcStretches(C); 
	DstrainEnergy();
	DDstrainEnergy();
	
	Matrix3d NN;
	
	NN = N1 * N1.transpose();
	VoigtFormat(NN,N1N1);
	
	NN = N2 * N2.transpose();
	VoigtFormat(NN,N2N2); 
	
	NN = N3 * N3.transpose();
	VoigtFormat(NN,N3N3); 
	
	NN = N1 * N2.transpose();
	VoigtFormat(NN,N1N2);
	
	NN = N1 * N3.transpose();
	VoigtFormat(NN,N1N3);
	
	NN = N2 * N1.transpose();
	VoigtFormat(NN,N2N1);
	
	NN = N2 * N3.transpose();
	VoigtFormat(NN,N2N3);
	
	NN = N3 * N1.transpose();
	VoigtFormat(NN,N3N1);
	
	NN = N3 * N2.transpose();
	VoigtFormat(NN,N3N2);
	
	if(!hypertype.compare("coupled"))
	{
		S1 = W1 / lambda1;
		S2 = W2 / lambda2;
		S3 = W3 / lambda3;
		
		Ct = (W11 - S1) / (lambda1 * lambda1) * N1N1 * N1N1.transpose() + W12 / (lambda1 * lambda2) * N1N1 * N2N2.transpose() + W13 / (lambda1 * lambda3) * N1N1 * N3N3.transpose() +
		     (W22 - S2) / (lambda2 * lambda2) * N2N2 * N2N2.transpose() + W12 / (lambda2 * lambda1) * N2N2 * N1N1.transpose() + W23 / (lambda2 * lambda3) * N2N2 * N3N3.transpose() +
		     (W33 - S3) / (lambda3 * lambda3) * N3N3 * N3N3.transpose() + W13 / (lambda3 * lambda1) * N3N3 * N1N1.transpose() + W23 / (lambda3 * lambda2) * N3N3 * N2N2.transpose();
		
		if(fabs(lambda1 - lambda2) > tol && fabs(lambda1 - lambda3) > tol && fabs(lambda2 - lambda3) > tol)
// 		if(lambda1 != lambda2 && lambda1 != lambda3 && lambda2 != lambda3 )
		{
			Ct += (S2 - S1) / (lambda2 * lambda2 - lambda1 * lambda1) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			      (S3 - S1) / (lambda3 * lambda3 - lambda1 * lambda1) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			      (S1 - S2) / (lambda1 * lambda1 - lambda2 * lambda2) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			      (S3 - S2) / (lambda3 * lambda3 - lambda2 * lambda2) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			      (S1 - S3) / (lambda1 * lambda1 - lambda3 * lambda3) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			      (S2 - S3) / (lambda2 * lambda2 - lambda3 * lambda3) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			       cout << "Entra en distintos" << endl;
		}
		else if ((fabs(lambda1 - lambda2) < tol) && (fabs(lambda2 - lambda3) < tol) )
// 		else if (lambda1 == lambda2 == lambda3) 
		{
			Ct += ((W22 - S2) / (2.0 * lambda2 * lambda2) - W12 / (2.0 * lambda1 * lambda2)) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			      ((W33 - S3) / (2.0 * lambda3 * lambda3) - W13 / (2.0 * lambda1 * lambda3)) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			      ((W11 - S1) / (2.0 * lambda1 * lambda1) - W12 / (2.0 * lambda1 * lambda2)) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			      ((W33 - S3) / (2.0 * lambda3 * lambda3) - W23 / (2.0 * lambda2 * lambda3)) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			      ((W11 - S1) / (2.0 * lambda1 * lambda1) - W13 / (2.0 * lambda1 * lambda3)) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			      ((W22 - S2) / (2.0 * lambda2 * lambda2) - W23 / (2.0 * lambda2 * lambda3)) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			      cout << "Entra en iguales" << endl;
		}
		else if (fabs(lambda1 - lambda2) < tol && fabs(lambda1 - lambda3) > tol)
// 		else if (lambda1 == lambda2 && lambda1 != lambda3)
		{
			Ct += ((W22 - S2) / (2.0 * lambda2 * lambda2) - W12 / (2.0 * lambda1 * lambda2)) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			      (S3 - S1) / (lambda3 * lambda3 - lambda1 * lambda1) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			      ((W11 - S1) / (2.0 * lambda1 * lambda1) - W12 / (2.0 * lambda1 * lambda2)) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			      (S3 - S2) / (lambda3 * lambda3 - lambda2 * lambda2) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			      (S1 - S3) / (lambda1 * lambda1 - lambda3 * lambda3) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			      (S2 - S3) / (lambda2 * lambda2 - lambda3 * lambda3) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			      cout << "Entra en l1=l2" << endl;
		}
		else if (fabs(lambda1 - lambda3) < tol && fabs(lambda1 - lambda2) > tol) 
// 		else if (lambda1 == lambda3 && lambda1 != lambda2 )
		{
			Ct += (S2 - S1) / (lambda2 * lambda2 - lambda1 * lambda1) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			      ((W33 - S3) / (2.0 * lambda3 * lambda3) - W13 / (2.0 * lambda1 * lambda3)) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			      (S1 - S2) / (lambda1 * lambda1 - lambda2 * lambda2) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			      (S3 - S2) / (lambda3 * lambda3 - lambda2 * lambda2) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			      ((W11 - S1) / (2.0 * lambda1 * lambda1) - W13 / (2.0 * lambda1 * lambda3)) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			      (S2 - S3) / (lambda2 * lambda2 - lambda3 * lambda3) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			      cout << "Entra en l1=l3" << endl;
		}
		else if (fabs(lambda2 - lambda3) < tol && fabs(lambda2 - lambda1) > tol)
// 		else if (lambda2 == lambda3 && lambda2 != lambda1)
		{
			Ct += (S2 - S1) / (lambda2 * lambda2 - lambda1 * lambda1) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			      (S3 - S1) / (lambda3 * lambda3 - lambda1 * lambda1) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			      (S1 - S2) / (lambda1 * lambda1 - lambda2 * lambda2) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			      ((W33 - S3) / (2.0 * lambda3 * lambda3) - W23 / (2.0 * lambda2 * lambda3)) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			      (S1 - S3) / (lambda1 * lambda1 - lambda3 * lambda3) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			      ((W22 - S2) / (2.0 * lambda2 * lambda2) - W23 / (2.0 * lambda2 * lambda3)) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			       cout << "Entra en l2=l3" << endl;
		}
// 		cout << "Ct= " << endl << Ct << endl << endl;
	}
	else //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	{
		VectorXd Ci(6);
		MatrixXd CioCi(6,6);
		Tensor4LSym IIc;
		
        calcCiTensor(C);
		IIc = selfOdot(CiTensor);
		VoigtFormat(IIc,CioCi);
		
		VoigtFormat(CSym.inverse(),Ci);
		
		Ctvol = J * (dWvol + J * ddWvol) * Ci * Ci.transpose() - 2.0 * J * dWvol * CioCi;
		
// 		Ctvol = J * (dWvol + J * ddWvol) * 
// 		        (N1N1 * N1N1.transpose() / pow(lambda1,4.0) + N1N1 * N2N2.transpose() / (lambda1 * lambda1 * lambda2 * lambda2) + N1N1 * N3N3.transpose() / (lambda1 * lambda1 * lambda3 * lambda3) +
// 		         N2N2 * N2N2.transpose() / pow(lambda2,4.0) + N2N2 * N1N1.transpose() / (lambda2 * lambda2 * lambda1 * lambda1) + N2N2 * N3N3.transpose() / (lambda2 * lambda2 * lambda3 * lambda3) +
// 		         N3N3 * N3N3.transpose() / pow(lambda3,4.0) + N3N3 * N1N1.transpose() / (lambda3 * lambda3 * lambda1 * lambda1) + N3N3 * N2N2.transpose() / (lambda3 * lambda3 * lambda2 * lambda2))
// 		        + J * dWvol * 
// 		        (N1N1 * N1N1.transpose() / pow(lambda1,4.0) + ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) / (lambda1 * lambda1 * lambda2 * lambda2) + ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) / (lambda1 * lambda1 * lambda3 * lambda3) +
// 		         N2N2 * N2N2.transpose() / pow(lambda2,4.0) + ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) / (lambda2 * lambda2 * lambda1 * lambda1) + ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) / (lambda2 * lambda2 * lambda3 * lambda3) +
// 		         N3N3 * N3N3.transpose() / pow(lambda3,4.0) + ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) / (lambda3 * lambda3 * lambda1 * lambda1) + ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) / (lambda3 * lambda3 * lambda2 * lambda2));
// 		
// 		cout << "Ctvol= " << endl << Ctvol << endl << endl;
		
// 		Siso11m = - Siso1 / lambda1m + (2.0 / 3.0 * (Wiso1 + lambda1m * Wiso11) - 1.0 / 3.0 * (lambda2m * Wiso12 + lambda3m * Wiso13)) / (lambda1 * lambda1);
// 		Siso22m = - Siso2 / lambda2m + (2.0 / 3.0 * (Wiso2 + lambda2m * Wiso22) - 1.0 / 3.0 * (lambda1m * Wiso12 + lambda3m * Wiso23)) / (lambda2 * lambda2);
// 		Siso33m = - Siso3 / lambda3m + (2.0 / 3.0 * (Wiso3 + lambda3m * Wiso33) - 1.0 / 3.0 * (lambda1m * Wiso13 + lambda2m * Wiso23)) / (lambda3 * lambda3);
// 		
// 		Siso12m = 2.0 * Siso1 / lambda2m + (2.0 / 3.0 * lambda1m * Wiso12 - 1.0 / 3.0 * (Wiso2 + lambda2m * Wiso22 + lambda3m * Wiso23)) / (lambda1 * lambda1);
// 		Siso13m = 2.0 * Siso1 / lambda3m + (2.0 / 3.0 * lambda1m * Wiso13 - 1.0 / 3.0 * (Wiso3 + lambda3m * Wiso33 + lambda2m * Wiso23)) / (lambda1 * lambda1);
// 		
// 		Siso21m = 2.0 * Siso2 / lambda1m + (2.0 / 3.0 * lambda2m * Wiso12 - 1.0 / 3.0 * (Wiso1 + lambda1m * Wiso11 + lambda3m * Wiso13)) / (lambda2 * lambda2);
// 		Siso23m = 2.0 * Siso2 / lambda3m + (2.0 / 3.0 * lambda2m * Wiso23 - 1.0 / 3.0 * (Wiso3 + lambda3m * Wiso33 + lambda1m * Wiso13)) / (lambda2 * lambda2);
//        
// 		Siso31m = 2.0 * Siso3 / lambda1m + (2.0 / 3.0 * lambda3m * Wiso13 - 1.0 / 3.0 * (Wiso1 + lambda1m * Wiso11 + lambda2m * Wiso12)) / (lambda3 * lambda3);
// 		Siso32m = 2.0 * Siso3 / lambda2m + (2.0 / 3.0 * lambda3m * Wiso23 - 1.0 / 3.0 * (Wiso2 + lambda2m * Wiso22 + lambda1m * Wiso12)) / (lambda3 * lambda3);
// 		
		////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
		/*
		Siso11 = (2.0 / 3.0 * lambda1m * Siso11m - 1.0 / 3.0 * (lambda2m * Siso12m + lambda3m * Siso13m)) / (lambda1 * lambda1);
		Siso22 = (2.0 / 3.0 * lambda2m * Siso22m - 1.0 / 3.0 * (lambda1m * Siso21m + lambda3m * Siso23m)) / (lambda2 * lambda2);
		Siso33 = (2.0 / 3.0 * lambda3m * Siso33m - 1.0 / 3.0 * (lambda1m * Siso31m + lambda2m * Siso32m)) / (lambda3 * lambda3);
		
		Siso12 = (2.0 / 3.0 * lambda2m * Siso12m - 1.0 / 3.0 * (lambda1m * Siso11m + lambda3m * Siso13m)) / (lambda2 * lambda2);
		Siso13 = (2.0 / 3.0 * lambda3m * Siso13m - 1.0 / 3.0 * (lambda1m * Siso11m + lambda2m * Siso12m)) / (lambda3 * lambda3);
		
		Siso21 = (2.0 / 3.0 * lambda1m * Siso21m - 1.0 / 3.0 * (lambda2m * Siso22m + lambda3m * Siso23m)) / (lambda1 * lambda1);
		Siso23 = (2.0 / 3.0 * lambda3m * Siso23m - 1.0 / 3.0 * (lambda2m * Siso22m + lambda1m * Siso21m)) / (lambda3 * lambda3);
		
		Siso31 = (2.0 / 3.0 * lambda1m * Siso31m - 1.0 / 3.0 * (lambda3m * Siso33m + lambda2m * Siso32m)) / (lambda1 * lambda1);
		Siso32 = (2.0 / 3.0 * lambda2m * Siso32m - 1.0 / 3.0 * (lambda3m * Siso33m + lambda1m * Siso31m)) / (lambda2 * lambda2);*/
		
// 		Siso11 = Jm / lambda1 * (2.0 / 3.0 * Siso11m - (lambda2m / lambda1m * Siso12m + lambda3m / lambda1m * Siso13m) / 3.0);
// 		Siso22 = Jm / lambda2 * (2.0 / 3.0 * Siso22m - (lambda1m / lambda2m * Siso21m + lambda3m / lambda2m * Siso23m) / 3.0);
// 		Siso33 = Jm / lambda3 * (2.0 / 3.0 * Siso33m - (lambda1m / lambda3m * Siso31m + lambda2m / lambda3m * Siso32m) / 3.0);
// 		
// 		Siso12 = Jm / lambda2 * (2.0 / 3.0 * Siso12m - (lambda1m / lambda2m * Siso11m + lambda3m / lambda2m * Siso13m) / 3.0);
// 		Siso13 = Jm / lambda3 * (2.0 / 3.0 * Siso13m - (lambda1m / lambda3m * Siso11m + lambda2m / lambda3m * Siso12m) / 3.0);
// 		
// 		Siso21 = Jm / lambda1 * (2.0 / 3.0 * Siso21m - (lambda2m / lambda1m * Siso22m + lambda3m / lambda1m * Siso23m) / 3.0);
// 		Siso23 = Jm / lambda3 * (2.0 / 3.0 * Siso23m - (lambda2m / lambda3m * Siso22m + lambda1m / lambda3m * Siso21m) / 3.0);
// 		
// 		Siso31 = Jm / lambda1 * (2.0 / 3.0 * Siso31m - (lambda3m / lambda1m * Siso33m + lambda2m / lambda1m * Siso32m) / 3.0);
// 		Siso32 = Jm / lambda2 * (2.0 / 3.0 * Siso32m - (lambda3m / lambda2m * Siso33m + lambda1m / lambda2m * Siso31m) / 3.0);
		
// 		cout << "Siso11-1= " << Siso11 << endl;
// 		cout << "Siso12-1= " << Siso12 << endl;
// 		cout << "Siso13-1= " << Siso13 << endl;
// 		cout << "Siso21-1= " << Siso21 << endl;
// 		cout << "Siso22-1= " << Siso22 << endl;
// 		cout << "Siso23-1= " << Siso23 << endl;
// 		cout << "Siso31-1= " << Siso31 << endl;
// 		cout << "Siso32-1= " << Siso32 << endl;
// 		cout << "Siso33-1= " << Siso33 << endl;
// 		cout << "---------------------"<< endl;
		
		Siso11 = 0.0;
		Siso22 = 0.0;
		Siso33 = 0.0;
		Siso12 = 0.0;
		Siso13 = 0.0;
		Siso21 = 0.0;
		Siso23 = 0.0;
		Siso31 = 0.0;
		Siso32 = 0.0;
		
		vector <double> Mu = getvMu();
		vector <double> Alpha = getvAlpha();
		
		for(int i=0;i<Mu.size();++i)
		{
			Siso11 += pow(lambda1,-4.0) * Mu[i] * ((4.0 / 9.0 * Alpha[i] - 4.0 / 3.0) * pow(lambda1m,Alpha[i]) + (2.0 / 3.0 + 1.0 / 9.0 * Alpha[i]) * (pow(lambda2m,Alpha[i]) + pow(lambda3m,Alpha[i])));
			Siso22 += pow(lambda2,-4.0) * Mu[i] * ((4.0 / 9.0 * Alpha[i] - 4.0 / 3.0) * pow(lambda2m,Alpha[i]) + (2.0 / 3.0 + 1.0 / 9.0 * Alpha[i]) * (pow(lambda1m,Alpha[i]) + pow(lambda3m,Alpha[i])));
			Siso33 += pow(lambda3,-4.0) * Mu[i] * ((4.0 / 9.0 * Alpha[i] - 4.0 / 3.0) * pow(lambda3m,Alpha[i]) + (2.0 / 3.0 + 1.0 / 9.0 * Alpha[i]) * (pow(lambda1m,Alpha[i]) + pow(lambda2m,Alpha[i])));
		}
		for(int j=0;j<Mu.size();++j)
		{
			Siso12 += pow(lambda1,-2.0) * pow(lambda2,-2.0) * Mu[j] * Alpha[j] * ( -2.0 / 9.0 * (pow(lambda1m,Alpha[j]) + pow(lambda2m,Alpha[j])) + 1.0 / 9.0 * pow(lambda3m,Alpha[j]) );
			Siso13 += pow(lambda1,-2.0) * pow(lambda3,-2.0) * Mu[j] * Alpha[j] * ( -2.0 / 9.0 * (pow(lambda1m,Alpha[j]) + pow(lambda3m,Alpha[j])) + 1.0 / 9.0 * pow(lambda2m,Alpha[j]) );
			
			Siso21 += pow(lambda2,-2.0) * pow(lambda1,-2.0) * Mu[j] * Alpha[j] * ( -2.0 / 9.0 * (pow(lambda2m,Alpha[j]) + pow(lambda1m,Alpha[j])) + 1.0 / 9.0 * pow(lambda3m,Alpha[j]) );
			Siso23 += pow(lambda2,-2.0) * pow(lambda3,-2.0) * Mu[j] * Alpha[j] * ( -2.0 / 9.0 * (pow(lambda2m,Alpha[j]) + pow(lambda3m,Alpha[j])) + 1.0 / 9.0 * pow(lambda1m,Alpha[j]) );
			
			Siso31 += pow(lambda3,-2.0) * pow(lambda1,-2.0) * Mu[j] * Alpha[j] * ( -2.0 / 9.0 * (pow(lambda3m,Alpha[j]) + pow(lambda1m,Alpha[j])) + 1.0 / 9.0 * pow(lambda2m,Alpha[j]) );
			Siso32 += pow(lambda3,-2.0) * pow(lambda2,-2.0) * Mu[j] * Alpha[j] * ( -2.0 / 9.0 * (pow(lambda3m,Alpha[j]) + pow(lambda2m,Alpha[j])) + 1.0 / 9.0 * pow(lambda1m,Alpha[j]) );
		}
// 		
// 		cout << "Siso11-2= " << Siso11 << endl;
// 		cout << "Siso12-2= " << Siso12 << endl;
// 		cout << "Siso13-2= " << Siso13 << endl;
// 		cout << "Siso21-2= " << Siso21 << endl;
// 		cout << "Siso22-2= " << Siso22 << endl;
// 		cout << "Siso23-2= " << Siso23 << endl;
// 		cout << "Siso31-2= " << Siso31 << endl;
// 		cout << "Siso32-2= " << Siso32 << endl;
// 		cout << "Siso33-2= " << Siso33 << endl;
// 		cout << "---------------------"<< endl;
		
		Ctiso = Siso11 * N1N1 * N1N1.transpose() + Siso12 * N1N1 * N2N2.transpose() + Siso13 * N1N1 * N3N3.transpose() +
		        Siso22 * N2N2 * N2N2.transpose() + Siso21 * N2N2 * N1N1.transpose() + Siso23 * N2N2 * N3N3.transpose() +
		        Siso33 * N3N3 * N3N3.transpose() + Siso31 * N3N3 * N1N1.transpose() + Siso32 * N3N3 * N2N2.transpose();
		
		if(fabs(lambda1 - lambda2) > tol && fabs(lambda1 - lambda3) > tol && fabs(lambda2 - lambda3) > tol)
// 		if(lambda1 != lambda2 && lambda1 != lambda3 && lambda2 != lambda3)
		{
			Ctiso += (Siso2 - Siso1) / (lambda2 * lambda2 - lambda1 * lambda1) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			         (Siso3 - Siso1) / (lambda3 * lambda3 - lambda1 * lambda1) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			         (Siso1 - Siso2) / (lambda1 * lambda1 - lambda2 * lambda2) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			         (Siso3 - Siso2) / (lambda3 * lambda3 - lambda2 * lambda2) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			         (Siso1 - Siso3) / (lambda1 * lambda1 - lambda3 * lambda3) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			         (Siso2 - Siso3) / (lambda2 * lambda2 - lambda3 * lambda3) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			       cout << "Entra en distintos" << endl;
		}
		else if ((fabs(lambda1 - lambda2) < tol) && (fabs(lambda2 - lambda3) < tol) )
// 		else if (lambda1 == lambda2 == lambda3) 
		{
			Ctiso += (Siso22 - Siso12) / 2.0 * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			         (Siso33 - Siso13) / 2.0 * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			         (Siso11 - Siso21) / 2.0 * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			         (Siso33 - Siso23) / 2.0 * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			         (Siso11 - Siso31) / 2.0 * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			         (Siso22 - Siso32) / 2.0 * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			      cout << "Entra en iguales" << endl;
		}
		else if (fabs(lambda1 - lambda2) < tol && fabs(lambda1 - lambda3) > tol)
// 		else if (lambda1 == lambda2 && lambda1 != lambda3)
		{
			Ctiso += (Siso22 - Siso12) / 2.0 * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			         (Siso3 - Siso1) / (lambda3 * lambda3 - lambda1 * lambda1) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			         (Siso11 - Siso21) / 2.0 * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			         (Siso3 - Siso2) / (lambda3 * lambda3 - lambda2 * lambda2) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			         (Siso1 - Siso3) / (lambda1 * lambda1 - lambda3 * lambda3) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			         (Siso2 - Siso3) / (lambda2 * lambda2 - lambda3 * lambda3) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			      cout << "Entra en l1=l2" << endl;
		}
		else if (fabs(lambda1 - lambda3) < tol && fabs(lambda1 - lambda2) > tol) 
// 		else if (lambda1 == lambda3 && lambda1 != lambda2 )
		{
			Ctiso += (Siso2 - Siso1) / (lambda2 * lambda2 - lambda1 * lambda1) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			         (Siso33 - Siso13) / 2.0 * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			         (Siso1 - Siso2) / (lambda1 * lambda1 - lambda2 * lambda2) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			         (Siso3 - Siso2) / (lambda3 * lambda3 - lambda2 * lambda2) * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			         (Siso11 - Siso31) / 2.0 * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			         (Siso2 - Siso3) / (lambda2 * lambda2 - lambda3 * lambda3) * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() ) ;
// 			      cout << "Entra en l1=l3" << endl;
		}
		else if (fabs(lambda2 - lambda3) < tol && fabs(lambda2 - lambda1) > tol)
// 		else if (lambda2 == lambda3 && lambda2 != lambda1)
		{
			Ctiso += (Siso2 - Siso1) / (lambda2 * lambda2 - lambda1 * lambda1) * ( N1N2 * N1N2.transpose() + N1N2 * N2N1.transpose() ) +
			         (Siso3 - Siso1) / (lambda3 * lambda3 - lambda1 * lambda1) * ( N1N3 * N1N3.transpose() + N1N3 * N3N1.transpose() ) +
			         (Siso1 - Siso2) / (lambda1 * lambda1 - lambda2 * lambda2) * ( N2N1 * N2N1.transpose() + N2N1 * N1N2.transpose() ) +
			         (Siso33 - Siso23) / 2.0 * ( N2N3 * N2N3.transpose() + N2N3 * N3N2.transpose() ) +
			         (Siso1 - Siso3) / (lambda1 * lambda1 - lambda3 * lambda3) * ( N3N1 * N3N1.transpose() + N3N1 * N1N3.transpose() ) +
			         (Siso22 - Siso32) / 2.0 * ( N3N2 * N3N2.transpose() + N3N2 * N2N3.transpose() );
// 			       cout << "Entra en l2=l3" << endl;
		}
// 		cout << "Ctiso= " << endl << Ctiso << endl << endl;
		
		Ct = Ctvol + Ctiso;
	}
// 		cout << "Ct= " << endl << Ct << endl << endl;
}

//***
